Breathing behavior involves a large network of neurons which are distributed throughout the nervous system in brain areas such as in the neocortex, cerebellum, amygdala, pons, medulla, and in the spinal cord.^[1,2] However, it is believed that one special region, being termed as the pre-Bötzinger complex (pre-BötC), seems to be essential for the generation and regulation of respiratory rhythms.^[3,4] Why can the breathing rhythm be exercised independently? It is because there are some inspiratory neurons acting like pacemaker cells in the pre-BötC.^[5] These clusters of neurons connect with each other, adjust breath and control inspiratory rhythm together.^[2,6,7]

Rhythm genesis in a neuronal system involves a quite complicated physical and chemical process, and at the same time includes various dynamic molecular mechanisms.^[1,2] Usually, neurons can connect with each other by chemical or electrical synapses to form a complex network to accomplish complicated dynamical behaviors.^[6–9] In order to explore the dynamical features of the neuronal firing activities, theories and methods in nonlinear dynamics have been explored in neuroscience to investigate the relationship between the neuronal behavior and dynamics.^[10–14] Among these methods, bifurcation analysis is an important one that has been heavily used in the research. Gu et al. discovered some types of bifurcation scenarios in neural firing rhythms in biological experiments.^[12,13] Duan et al. investigated the regular bursting generation and pattern transitions in the single or two-cell model network, and found that the pattern transitions are closely related to the codimension-one and -two bifurcations.^[15–17] Also, codimension-2 Bautin bifurcation of synchronous solution of a coupled FHN neural system with delay was investigated by Zhen.^[18] Abundant discharge of neurons activity, from the view point of dynamics, is caused by various dynamic bifurcations.^[14,19–21]

A novel bursting mode has been found experimentally in the pre-Bötzinger neurons, which is referred to as mixed bursting.^[22,23] Dunmyre suggested that such bursting may result from the interactions of the persistent sodium and calcium-activated nonspecific cationic currents (NaP and CAN currents).^[24] Different from the mixed-mode oscillation,^[25] mixed bursting is a special periodic oscillation solution of multi-timescale systems characterized by containing different types of short bursts in one period.^[26] Mixed bursting solutions observed in biological models or experiments contain more dynamical information and are worthy of further study. The generation mechanism of mixed bursting is very complex which involves disparate time scales. Using fast-slow decomposition,^[27] Wang and Rubin^[26] showed that mixed busting can exist in two-timescale systems. Rubin et al.^[28] presented an approach that can guide model development and tuning to achieve desired qualitative and quantitative solution properties.

Voltage-gated potassium current is important and it can trigger endogenous oscillations in neurons. Shevtsova et al.^[29] showed that increased potassium concentration could increase the bursting frequency and decrease the bursting duration, and at a higher potassium concentration, bursting switches to tonic firing. The delayed-rectifier potassium could trigger rhythmic bursting activity at single pacemaker neurons and network population.^[30–32] However, although the effects of potassium currents on the bursting and spiking were demonstrated, little attention has been paid to its effects on the mixed bursting. The focus of this paper is on the analysis of the dynamical mechanisms underlying such mixed bursting with the varying of the potassium current.

This paper is organized as follows. In Section 2, we introduce the model of coupled pre-BötC neurons.^[6,7] In Section 3, based on the model introduced in Section 2, mixed bursting patterns and their transitions are studied. We explore the influences of the potassium concentration on the mixed bursting by fast-slow bifurcation decomposition. The conclusion is given in Section 4.

The model of coupled excitable pre-BötC inspiratory neurons introduced by Butera^[6,7] is described as follows:

When the whole cell capacitance is changed, ε/τ_h(V_i) is far less than 1/τ_n(V_i). The derivative of h_i is close to 0, which means that the changing rate of h_i is much lower than that of other variables. Therefore, h₁ and h₂ can be regarded as slow variables and the model can be considered as a fast-slow decomposition with Eqs. (1), (3), and (4) composing the fast subsystem, and Eq. (2) being the slow subsystem. In what follows, we set ε = 6, and the coefficient of derivative of h_i, . Then h₁ and h₂ are approximately equal, and we set h₁ ≈ h₂ = h. The bifurcation is performed by software XPPAUT.^[33]

Bursting is one of the most important firing patterns responsible for the information transitions in the pre-Bötzinger complex. In the two-coupled cells network, bursting oscillations can be divided into two types, namely, the in-phase oscillations and the anti-phase oscillations.^[19] We obtain the in-phase oscillation when the initial values for neuron 1 and neuron 2 are the same. Otherwise, we obtain the anti-phase oscillations. The bursting pattern depends on different ion currents, especially on the ion conduction, for example g_K. In the following, considering g_K as a changing parameter, we will study the mixed anti-phase burstings, especially their generation and transitions through bifurcation analysis.

The coupled neurons are not exactly synchronous but can exhibit the similar bursting pattern under the same set of parameters. So we focus on V₁ to illustrate our analysis and results. The bursting oscillations of V₁ corresponding to different values of g_K are shown in Fig. 1. The neuron activity undergoes a transition through silence, regular bursting, mixed bursting, and regular bursting finally with the increase of parameter g_K. At the same time, the frequency of the bursting oscillation is increased.

Fig. 1.

	Figure Option ViewDownloadNew Window
	Fig. 1. Time series of system (1)–(4) corresponding to different values of gK: (a) gK = 4.4 nS, (b) gK = 5.2 nS, (c) gK = 5.92 nS, (d) gK = 6.6 nS, (e) gK = 7.1 nS, (f) gK = 7.5 nS.

3.1. Bifurcation analysis for silence

When g_K = 4.4 nS, the coupled neurons are in a silent state, as shown in Fig. 1(a). Figure 2(a) shows the trajectory of system (1)–(4) on the (n₁, V₁) plane from which we can see that the trajectory starts from point P₀ and ends at point P₁. Bifurcation analysis of the fast subsystem is shown in Fig. 2(b) with the slow variable h = h₁ = h₂ being a bifurcation parameter. The three S-shaped curves are composed of the equilibrium points with fold bifurcation point on the knees (F₁ and F₂), two subcritical Hopf bifurcations (Subh₁ and Subh₂) on the upper branch where the unstable focus transits to a stable focus. Unstable limit cycles emanate from the subcritical Hopf bifurcation points Subh₁ and Subh₂. LPC₁ and LPC₂ refer to the fold bifurcations of limit cycles where stable and unstable limit cycles coincide and disappear. The nullcline of variable h (h = h₁) is also appended. The h-nullcline interacts with the upper branch of the S-shaped curve at the point of stable focus P₁, right to the subcritical Hopf bifurcation points. So the trajectory of system (1)–(4) starting from point P₀ heads to the lower stable state and then transits to the upper state by the fold bifurcation point F₁ on the anticlockwise. Then the trajectory dwells around the upper branch of the stable focus and heads to the fixed point P₁ lastly. So the silence potential is about −21.97 mV, as shown in Fig. 1.

Fig. 2.

Figure Option ViewDownloadNew Window

Fig. 2. (color online) (a) Trajectory of system (1)–(4) on the (n1,V1) plane for gK = 4.4 nS. The olive curve is the trajectory of the whole system, and the black and pink curves represent the V1- and n1-nullcline. (b) Bifurcation analysis of the fast subsystem with the slow variable as a parameter. The red solid curve and dashed curve describe stable and unstable limit cycles, respectively, the olive curve shows the trajectory of the whole system, and the blue curve is the h1-nullcline.

3.2. Bifurcation analysis for regular bursting

As shown in Fig. 1(b), the silence state of the two-coupled cells changes to a firing state when g_K increases to 5.2 nS. The bifurcation structure of the fast subsystem is projected onto the (h, V₁) plane, as shown in Fig. 3(a). The meanings of labels for points and curves are the same as those shown in Fig. 2(b). Note that two families of periodic orbits emanate from the distinct subcritical Hopf bifurcations labeled Subh₁ and Subh₂ on the upper branch of the S-shaped curve. The red solid and dashed curves represent the maximum and minimum values of stable and unstable limit cycles, respectively. The h-nullcline intersects the upper branch of the S-shaped curve at point P, which is an unstable focus that lies on the left beside the two subcritical Hopf bifurcations, as shown in Fig. 3(b).

The bistability is formed by the stable nodes on the down state and the stable focuses on the up state. The trajectory of the system transits from the down rest state to the upper steady state by fold bifurcation (F₁). Then the trajectory oscillates around the stable focuses with damped amplitude until it passes through the Hopf bifurcation as a result of the slow passage effect. After that the amplitude of the trajectory increases gradually due to the repellent of the unstable focus. Finally, the trajectory transits to the lower stable state and completes a periodic oscillation. The two-coupled neurons exhibit the subHopf/subHopf type bursting.^[20]

Fig. 3.

	Figure Option ViewDownloadNew Window
	Fig. 3. (color online) (a) Fast-slow decomposition analysis with gK = 5.2 nS. (b) Enlargement of (a).

3.3. Bifurcation analysis for mixed bursting

As g_K increases to 5.92 nS, the system exhibits mixed bursting, which contains four different types of short bursts in one period, as shown in Fig. 1(c). Bifurcation analysis on the short bursts is shown in Fig. 4. Among them, the active phase of the first bursting starts and ends with the subcritical Hopf bifurcation. The bistability is formed by the stable nodes on the down state and the stable limit cycles on the up state. The rest state of the bursting transits to the up state by fold bifurcation F₁, and the firing state transits to the down state by fold bifurcation of limit cycle LPC₁. So the system exhibits the fold/fold limit cycle bursting, as shown in Fig. 4(a).

The other two bursting patterns (shown in Figs. 4(b) and 4(c)) are the subHopf/subHopf type bursting according to the classification scheme of Izhikevich,^[14] which is similar to that in Fig. 3. There is no attraction of limit cycles, and the bistability is still formed by the stable nodes and the stable focuses. For the fourth pattern, the main difference is that there are many small oscillations when the trajectory passes through the slow passage formed by the subcritical Hopf bifurcation. The bistability is formed by the stable nodes and the stable limit cycles, but the attraction of the stable limit cycles becomes weak due to the delay effect of the subHopf bifurcation. The bursting is the subHopf/fold limit cycle type via fold/fold limit cycle hysteresis loop, as shown in Fig. 4(d). The subHopf/fold limit cycle bursting with small oscillation is a transit state between the subHopf/subHopf bursting (Figs. 4(b) and 4(c)) and the fold/fold limit cycle bursting with big oscillation (Fig. 4(a)). Thus there are three different types of short bursts in the mixed bursting.

When g_K further increases to 6.6 nS, two types of short bursts occur alternatively, as shown in Fig. 1(d). Bifurcation analysis of these two short burstings is shown in Fig. 5. For the first one, the rest state transits to a firing state through fold bifurcation F₁, and then the firing state transits to the rest state by fold bifurcation of limit cycles LPC₁, as shown in Fig. 5(a). Thus it is the fold/fold limit cycle bursting (the same as that shown in Fig. 4(a)). According to Fig. 5(b), the second bursting is subHopf/fold limit cycle bursting via fold/fold limit cycle hysteresis loop, which is the same as that shown in Fig. 4(d).

If g_K increases to 7.1 nS, there are many different bursting patterns and all of them are irregular (Fig. 1(e)). Some of them have the same characteristics with the fold/fold limit cycle bursting, as shown in Figs. 6(a), 6(b), 6(d), and 6(e). The others have the same characteristics with the subHopf/fold limit cycle bursting via fold/fold limit cycle hysteresis loop, as shown in Fig. 6(c).

Fig. 4.

	Figure Option ViewDownloadNew Window
	Fig. 4. (color online) Fast-slow decomposition analysis of the fast subsystem for mixed bursting with gK = 5.92 nS: (a) fold/fold limit cycle bursting, (b) and (c) subHopf/subHopf type bursting, (d) subHopf/fold limit cycle bursting. The meanings of labels and curves are the same as those in Fig. 3.

Fig. 5.

	Figure Option ViewDownloadNew Window
	Fig. 5. (color online) Fast-slow decomposition analysis for mixed bursting with gK = 6.6 nS: (a) fold/fold limit cycle bursting, (b) subHopf/fold limit cycle bursting. The meanings of labels and curves are the same as those in Fig. 3.

Fig. 6.

	Figure Option ViewDownloadNew Window
	Fig. 6. (color online) Fast-slow decomposition analysis for mixed bursting with gK = 7.1 nS: (a), (b), (d), and (e) fold/fold limit cycle bursting, (c) subHopf/fold limit cycle bursting. The meanings of labels and curves are the same as those in Fig. 3.

3.4. Two-parameter bifurcation analysis

When g_K = 7.5 nS, the system transits to a regular bursting, as shown in Fig. 1(f). The bifurcation analysis shown in Fig. 7 implies that the system exhibits the fold/fold limit cycle bursting.

Fig. 7.

	Figure Option ViewDownloadNew Window
	Fig. 7. (color online) Fast-slow decomposition analysis for regular bursting with gK = 7.5 nS. The meanings of labels and curves are the same as those in Fig. 3.

Mixed bursting is a kind of irregular bursting involving two-timescale phenomenon.^[26] Many factors, such as the types of bifurcation and relative positions of bifurcation points, may have an influence on the generation of such bursting. We will take advantage of two-parameter bifurcation in our analysis to identify how variations in g_K affect the fast variable dynamics and further lead to the existence of complex bursting solutions. When g_K increases from 5.2 nS to 7.5 nS, the Hopf bifurcation moves from left to right but the fold bifurcation (F₁) is almost at the same position, as shown in Fig. 8(a). As we have seen, mixed bursting contain different types of short bursts that can be observed during the variation of g_K. Considering h and g_K as parameters, we show the two-parameter bifurcation of the fast subsystem in Fig. 8(b). Changes of the relative position for the Hopf bifurcation curves (subh, related to Subh₁ and Subh₂, two subHopf bifurcation curves overlapped), the fold bifurcation curve of equilibrium points (f, related to right fold bifurcation F₁), and the fold bifurcation of limit cycles (l, related to LPC₁ and LPC₂) are shown in Fig. 8(b). The relative positions of the bifurcation curves have changed in region II, which is related to the mixed bursting (see Fig. 8(b)). When the parameter g_K lies in region I, the coupled neurons keep in the silence state. When g_K lies in region II, regular subHopf/subHopf bursting appears. When g_K lies in region III, different bursting types appear simultaneously, which give rise to the mixed bursting. When the parameter g_K lies in region VI, mixed bursting disappears and the system exhibits regular bursting again. See Fig. 8(b) and Table 1 for details.

Table 1.

Table 1.

Table 1. Approximate parameter regions for different neuronal activities. . Region Parameter range Neuronal activities region I 0 < gK ≤ 4.4 nS silence state region II 4.4 nS < gK < 5.87 nS regular subHopf/subHopf bursting region III 5.87 nS ≤ gK ≤ 7.45 nS mixed bursting region VI 7.45 nS < gK < 20 nS regular fold/fold limit cycle bursting

Table 1. Approximate parameter regions for different neuronal activities. .

Fig. 8.

Figure Option ViewDownloadNew Window

Fig. 8. (color online) (a) Fast-slow decomposition analysis of the fast subsystem with gK = 5.2 nS (red), 5.92 nS (blue), and 7.5 nS (green) respectively. (b) Two-parameter bifurcation with parameters gK and h. The black curve represents the Hopf bifurcation curve, the red curve represents the fold bifurcation of equilibrium, and the blue curve represents the fold bifurcation of limit cycles. CP is the codimension-2 cusp bifurcation.

Using both fast-slow decomposition and two-parameter bifurcation analysis, we investigate the generation mechanisms of the mixed bursting. We show that the pattern of mixed bursting in two-coupled pre-Bötzinger complex cells varies with the increase of g_K. Mixed bursting is widespread and complex. Dunmyre^[24] found in a self-coupled pre-Bötzinger the mixed patterns of square-wave and DB bursts (bursts featuring depolarization block), which matched well with the activity that was observed in the experimental preparations.^[23] In this paper, we find the mixed bursting composed of the regular fold/fold limit cycle bursting, subHopf/subHopf bursting, or subHopf/fold limit cycle bursting in the coupled pre-Bötzinger complex model.

Wang et al.^[26] studied the mixed bursting of a two-compartment model of a pre-Bötzinger complex neuron featuring a persistent sodium current I_NaP, a nonspecific cation CAN current, as well as the intracellular calcium oscillations that modulate the CAN current. Their investigations showed that the appearance of mixed bursting in pre-Bötzinger neurons depends heavily on the contribution of the dendritic calcium oscillations. We investigate the mixed bursting raised by the change of the potassium current (parameter g_K) in the coupled pre-Bötzinger complex model, which involves persistent sodium current I_NaP, but has no nonspecific cation CAN currents. Our results show that the nonspecific cation CAN current is not necessary for the generation of mixed bursting. The results also indicate that the mixed bursting is a heterogeneous transition state. Our analysis of mixed bursting patterns can be extended and applied to investigate other rhythmic neuronal systems with a two-timescale forcing term.